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RESUME 


On  derive  un  systeme  ferme  de  trois  equations  simultanees  aux 
derivees  partielles  pour  resoudre  l'intensite  moyenne  et  la  variance 
de  l'intensite  de  faisceaux  laser  focalises  se  propageant  dans  un  milieu 
turbulent.  Ces  equations  sont  uniformement  valides  a  tous  les  niveaux 
de  scintillation.  Des  exemples  de  solutions  calculees  numeriquement 
sont  en  excellent  accord  avec  des  mesures  obtenues  dans  1' atmosphere 
et  en  laboratoire.  On  derive  de  ces  equations  une  expression  analytique 
de  la  diffusion  par  la  turbulence  et  on  demontre  que  celle-ci  est  plus 
precise  que  la  formule  generalement  employee  et  tiree  du  concept  du 
diametre  de  coherence  p0.  Dans  la  limite  des  scintillations  faibles, 
la  solution  pour  la  variance  de  l'intensite  recoupe  les  resultats 
classiques  de  la  theorie  de  perturbation.  (NC) 


ABSTRACT 


A  closed  set  of  three  simultaneous  partial  differential  equations 
is  derived  for  the  solution  of  the  average  irradiance  and  the  irradiance 
variance  of  focused  laser  beams  in  a  turbulent  medium.  The  equations  are 
uniformly  valid  for  arbitrary  scintillation  levels.  Examples  of  solutions 
calculated  by  finite  difference  techniques  agree  very  well  with  experi¬ 
mental  data  in  the  atmosphere  and  in  the  laboratory.  An  analytic  expres¬ 
sion  for  the  turbulence-induced  beam  spreading  is  obtained  and  shown  to 
be  more  precise  than  the  commonly  used  formula  based  on  the  coherence 
diameter  p0.  In  the  weak  scintillation  limit,  the  irradiance  variance 
solution  agrees  with  the  classical  perturbation  results.  (U) 
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1.0  INTRODUCTION 

Atmospheric  turbulence  can  seriously  affect  the  propagation  of 
a  laser  beam.  The  refractive  index  turbulence  induces  phase  and  ampli¬ 
tude  fluctuations  which  cause  the  beam  to  scintillate  or  break  up  in 
several  random  patches,  to  lose  its  spatial  coherence,  to  wander  about 
its  axis,  and  to  spread  out.  Numerous  analytical  treatments  of  these 
phenomena  have  already  been  published.  The  analytical  Rytov  method, 
fully  documented  in  Refs.  1  and  2,  gives  excellent  results  in  the  weak 
scintillation  regime,  i.e.  at  small  propagation  distance  and/or  turbu¬ 
lence  strength.  However,  it  is  well-known  that  the  Rytov  method  fails 
when  irradiance  fluctuations  reach  a  certain  level:  theory  predicts  an 
unlimited  increase  of  the  irradiance  variance  whereas  experiments 
definitely  show  saturation.  Many  recent  models  deal  with  this  problem 
and  predict,  with  variable  degrees  of  accuracy,  the  phenomenon  of  satu¬ 
ration.  A  first  approach,  known  as  the  renormalization  technique,  is 
described  in  Ref.  2.  For  example,  it  is  used  by  De  Wolf  (Refs.  3-6) 
to  predict  the  saturation  of  the  irradiance  variance.  A  second  approach 
consists  in  solving  the  differential  equation  for  the  fourth-order 
statistical  moment  of  the  complex  electric  field  closed  by  the  local 
application  of  the  method  of  small  perturbations.  This  method  is 
reviewed  in  Ref.  7  and  examples  of  solutions  showing  saturation  are 
given  in  Refs.  7-11.  A  third  approach  is  based  on  the  extension  of 
the  Huygens-Fresnel  principle  to  turbulent  media  (Ref.  12).  Saturation 
results  derived  from  various  applications  of  this  principle  may  be  found 
in  Refs.  13-17. 

The  models  listed  in  the  preceding  paragraph  are  generally  well 
corroborated  by  measurements.  However,  they  all  require  lengthy 
calculations  to  compute  the  average  irradiance  and  the  irradiance  variance 
profiles  of  finite  laser  beams.  One  has  either  to  solve  a  fourfold 
partial  differential  equation  or  to  evaluate  multiple  integrals  at  each 
spatial  point.  In  particular,  these  models  are  not  suitable  for  the 
calculation  of  nonlinear  propagation  in  the  presence  of  turbulence. 
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such  as  thermal  blooming.  To  fill  this  need,  the  present  report  pro¬ 
poses  a  calculation  method  of  beam  propagation  in  turbulence  which  is 
numerically  simple,  fast  and  uniformly  applicable  to  arbitrary  scin¬ 
tillation  levels. 

In  Ref.  18,  we  developed  a  model  in  the  form  of  linear,  second- 
order  partial  differential  equations  for  the  first-  and  second-order 
moments  of  the  complex  wave  amplitude.  These  equations  can  be  easily 
solved  numerically  to  calculate  the  average  irradiance  and  the  irradiance 
variance  profiles  of  laser  beams  in  turbulence.  The  approach  is  similar 
to  that  used  in  the  calculation  of  turbulent  shear  flows  (see  for 
instance  Chap.  5  of  Ref.  19).  The  differential  equations  for  the  statis¬ 
tical  moments  of  the  complex  wave  amplitude  are  derived  from  the  sto¬ 
chastic  wave  equation.  Since  the  latter  is  quasi-linear  owing  to  the 
nonhomogeneous  refractive  index,  the  resulting  system  of  equations  is 
unclosed.  To  resolve  this  difficulty,  expressions  are  sought  to  relate 
the  unknown  or  higher  order  moments  to  the  first-  and  second-order 
complex -amplitude  moments.  In  Ref.  18,  these  closure  relations  were 
derived  from  mostly  empirical  considerations  but  the  accuracy  of  the 
resulting  solutions  was  so  consistent  that  it  prompted  us  to  investigate 
further  into  this  approach. 

The  closure  relations  necessary  to  solve  the  average  irradiance 
were  rederived  in  Ref.  20  in  a  fully  consistent  manner  from  the  governing 
wave  equation.  However,  one  empirically  based  relation  is  still  needed 
to  predict  the  higher  order  irradiance  variance.  The  empirical  input 
appears  in  the  form  of  one  universal  complex  constant.  To  eliminate 
this  constant,  we  would  have  to  consider  the  equation  for  the  covariance 
of  the  complex  amplitude  and,  thus,  go  back  to  the  complexity  of  at 
least  fourfold  partial  differential  equations  which  we  set  out  to  avoid. 
Hence,  we  believe  that  the  level  of  empiricism  remaining  in  the  model  is 
acceptable  and  certainly  well  compensated  for  by  the  simplicity  it 
affords . 
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The  results  of  Ref.  20  are  strictly  applicable  to  collimated 
beams.  In  this  report,  we  extend  the  model  to  include  focalization. 
Section  2.0  is  a  short  review  of  the  theoretical  background.  Section 
3.0  gives  the  closure  relations.  Section  4.0  compares  solutions  with 
data  taken  in  a  simulation  experiment  and  in  the  atmosphere.  Finally, 
section  5.0  presents  simplified  analytical  solutions  valid  under  practical 
operating  conditions. 

This  work  was  performed  at  DREV  between  October  78  and  November 
79  under  PCN  33B07,  Atmospheric  Propagation  of  Laser  Beams. 


2.0  THEORETICAL  BACKGROUND 

Our  mathematical  model  for  optical  wave  propagation  is  described 
in  Refs.  18  and  20.  In  short,  the  scalar  electric  field  E  of  a  monochroma¬ 
tic  wave  propagating  in  the  z  direction  under  negligible  polarization 
effects  satisfies  the  equation 


V  E 


k2N2 


E 


[1] 


where  £  is  the  three-dimensional  gradient  operator,  k  =  n0w/c  is  the 
optical  wave  number,  n0  is  the  unperturbed  index  of  refraction,  w  is 
the  optical  angular  frequency  of  the  source,  c  is  the  speed  of  light 
in  free  space,  and  N  is  the  instantaneous  random  index  of  refraction. 
The  electric  field  E  is  written  as  follows: 


i 


E  -  A  exp[ik(4>+z)  -  imt]  , 


[2] 
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where  A  and  4  will  be  specified  below.  Upon  substituting  eq .  2  for  E 
in  eq.  1  and  making  the  paraxial  approximation  valid  for  propagation  in 
atmospheric  turbulence,  we  obtain  after  separation 


[3] 


[4] 


where  is  the  gradient  operator  in  the  plane  normal  to  the  0-z  axis, 
and 


V  =  £<j>  . 


[5] 


In  obtaining  eq.  3,  it  was  assumed  that  (N-n0)/n0  <<1  which  is  well 
verified  in  atmospheric  turbulence. 

The  separation  into  eqs.  3  and  4  is  different  from  common  practice. 
Consequently,  functions  A  and  <t>  differ  from  classical  amplitude  and 
phase.  The  advantage  of  this  particular  separation  is  that  the  equation 
for  V  is  independent  of  A.  It  is  easy  to  verify  that  eq.  3  is  the 
paraxial  form  of  the  eikonal  equation  of  geometrical  optics.  Hence, 
the  surfaces  41  =  constants  are  the  geometrical-wave  fronts  and  Vfz.r) 
is  the  component,  in  the  transverse  plane,  of  the  unit  vector  parallel 
to  the  geometrical  ray  passing  through  the  point  (z,r_).  V(z,r)  can 
also  be  interpreted  as  the  vector  angle  subtended  by  the  ray  at  that 
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point  and  will  be  referred  to  as  the  angl e-of -arrival .  -  mully,  from 

eq.  2,  it  follows  that  A  is  the  complex  amplitude  defined  on  the  geo¬ 
metrical  phase  fronts.  This  complex  amplitude  embodies  the  diffractional 
phase . 


Since  the  refractive  index  N  is  a  random  function,  eqs.  3  and  4  are 
stochastic  equations.  There  is  no  known  general  method  of  solving  the 
random  functions  and  A.  Even  if  that  could  be  achieved,  it  would  yield 
much  more  information  than  required  in  practice.  Here,  the  stochastic 
eqs.  3  and  4  are  used  to  obtain  the  differential  equations  for  the 
statistical  moments  of  V  and  A.  We  limit  ourselves  to  the  first-  and 
second-order  moments  which  lend  themselves  to  straightforward  physical 
interpretations. 

The  random  functions  are  written  as  sums  of  an  average  and  a 
fluctuating  part,  i.e.: 


N  =  <N>  +  n;  <n>  =  0, 


[6] 


V  =  <V>  +  v:  <v>  =  0, 

/V  7  /«*»«  7 


[7] 


A  =  <A>  +  a;  <a>  =  0, 


[8] 


where  the  angular  brackets  denote  ensemble  averaging.  The  equations 
for  the  moments  are  obtained  by  taking  the  ensemble  average  of  the 
stochastic  equations  for  _V,  A,  AA  and  AA*  derived  from  the  governing 
eqs.  3  and  4.  We  find: 


m 


■  <*' 
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(  —  +  <V>  •  V  \  <v> 

\3z  / 


=  -  1/2  ^<<v*x>  +  X(<N>  -  n0)/nc 


[9] 


/—  +  <V>*V  j  <A>  +  1/2  <A>V*<V>  -  — —  V2<A>  =  -V*<ya>+  l/2-  aV-^>) 
\3z  /  2k 

[10] 


f  —  +  <V>  -V  \  <aa>  + 
V z  ~ 


i  i 

'aa>V*<V>  -  —  7‘<aa>  +  —  <Va*Va> 


2k  k 

=  -  V*<vaa>  -  2<av>*V<A>  -  <A><aV»v> 


[11] 


(—  +  <V>*V  \<aa*>  + 
3z  J 


<aa*>V*<V>  -  — —  [<a*V2a>  -  <aV2a*>] 
~~  2k 


-  V «<vaa*>  -  <av> • V<A>*  -  <a*v>*V<A> 


-  1/2  <A>*<aV,*v>  -  1/2  <A><a*X_* v>  ( 


[12] 


where  the  superscript  *  denotes  a  complex  conjugate. 

The  system  of  eqs.  10-12  for  <A>,  <aa>  and  <aa*>  is  mathemati¬ 
cally  unclosed  insofar  as  it  contains  more  unknowns  than  equations, 
namely  the  moments  <j/a>,  <vna>,  <vaa*>  and  <aV,-„v>.  In  principle,  one 
could  derive  equations  for  these  moments  using  a  procedure  analogous 
to  that  used  in  obtaining  eqs.  10-12.  However,  one  would  quickly 
realize  that  this  technique  leads  to  equations  involving  still  higher 
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order  and  unknown  moments.  This  is  the  classical  closure  problem  which 
is  always  encountered  in  the  treatment  of  statistical  phenomena  governed 
by  nonlinear  and/or  quasi-linear  stochastic  equations  such  as  eqs.  3  and  4. 
Of  course,  there  is  no  exact  method  of  solution  for  this  problem  since 
the  complete  mathematical  model  contains  an  infinite  number  of  equations. 
Hence,  workable  models  require  closing  the  hierarchy  of  equations  after 
moments  of  a  given  order.  This  can  be  accomplished  only  through 
approximations  regarding  the  higher  order  moments  for  which  the  equations 
have  been  left  out.  In  this  paper,  we  derive  four  relations  to  express 
cya>,  <vaa> ,  <vaa*>  and  <a,V*y>  as  functions  of  <A> ,  <aa>  and  <aa*>. 

3.0  CLOSURE  RELATIONS 

The  derivation  of  the  closure  relations  is  based  on  three  prin¬ 
cipal  hypotheses. 

I.  The  fluctuating  angle-of-arrival  y  is  statistically  homoge¬ 
neous  and  isotropic  in  the  plane  transverse  to  the  0-z  axis.  This 
approximation  follows  from  the  hypothesis  that  the  refractive  index  is 
statistically  homogeneous  and  isotropic  and  from  the  paraxial  approxima¬ 
tion.  Indeed,  since  the  medium  is  homogeneous  and  isotropic  and  since 
the  average  ray  angle  <V>  is  small,  the  rays  reaching  a  plane  z  traverse 
statistically  equivalent  paths.  Hence,  the  covariance  function  of  v 
should  depend  only  on  the  relative  positions  of  the  observation  points 
in  plane  z. 

II.  The  random  complex  amplitude  and  angle-of-arrival  are  only 
weakly  correlated.  This  approximation  is  based  on  the  observation  that 
the  fluctuating  complex  amplitude  a(z,y)  is  the  result  of  repeated 
interactions  with  the  angle-of-arrival  over  the  complete  propagation 
path  between  0  and  z.  Hence,  a(z,r)  depends  on  many  processes  indepen¬ 
dent  of  the  local  angle-of-arrival  v,(z,r).  Moreover,  eq.  3  shows  that 


.■nmuM 
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v  is  independent  of  "a".  Therefore,  the  hypothesis  of  weak  correlation 
appears  very  plausible.  More  specifically,  we  assume 


< v  a  >  << 

2 


<V  •  v: 


1/2 


[13a] 


r  1/2 

<v  v  a  >  «  <v*v>  <I> 

~1~2  2  ~  ~ 


[13b] 


<v  v  a  a  >  -  <v  v  ><a  a  >  +  (|terms|<<  <v*v><I>}  [13c] 

~1~2  1  2  ~1~2  12 


etc .  .  .  , 


where  <I>  =  <A><A>*  +  <aa*>  is  the  average  irradiance.  The  subscripts 
1  and  2  refer  to  the  separation  points  of  the  covariance  functions; 
where  no  subscripts  are  used,  the  midpoint  between  1  and  2  is  taken. 

III.  The  covariance  functions  involving  the  fluctuating  complex 
amplitude  are  quasi -homogeneous  and  quasi-isotropic  in  the  transverse 
plane.  More  specifically,  we  set 


<a(z  ,r  )a(z  ,r  )>=  F  (z  ,z  ; 
112  2  aa1-  i  2 


r  +r 

Zl 


3aa 


[z,  >z-,\r-r  |)  ,  [14a] 


'1  '-h2 
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<a(z  ,  r  )a*(z  ,r  )> 
1  ~1  2  ~2 


r  +£ 

f  .*(**;  ~1  ~2  )  g„.(z 


aax  i  2 


’aa 


[  14b  J 


Under  hypothesis  I,  eq.  9  becomes 


—  <J/>  +  <)£>•&<£*  =  /V(<N>-n0)/n0,  [15] 

3z 

which  shows  that  the  average  ray  angle  <V>  is  independent  of  the  turbu¬ 
lent  fluctuations.  The  solution  of  eq.  15  for  an  initially  spherical 
phase  front  with  radius  of  curvature  F  in  a  statistically  homogeneous 
medium,  i.e.  <N>  =  constant,  is  given  by 


<V>  =  r/(z-F). 


[16] 


The  equation  for  the  fluctuating  complex  amplitude  is  obtained 
by  subtracting  eq.  10  from  eq.  4.  We  find,  using  eq.  16 


x* 

—  a  +  — —  *Va  +  — - -  -  —  V2a  =  g(z,r),  [17] 

3z  (z-F)  "  (z-F)  2k 


where 


<A>V «v 

g(z,r)  =  -v* V/<A>  -  - —  -  V- [x,a-<xa>]  +  l/2[aV*v-<a^-v>]  [18] 
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An  implicit  solution  to  eq.  17  with  the  boundary  condition 


a (0 ,r ) 


lim 
|r  I-*-00 


a(z,r)  =  0 


[19] 


is  given  by 


a(z,r)  =  L(z,x;  u.s,)  g(u,s,)  , 


[20] 


where  L(z,r;  u.s^)  is  an  integral  operator  defined  as  follows:’*" 


L(z»r;uf^)g(u,^) 


—  [  — —  [I  d^s  .(uoJexp 
2*)o  ^  U 


ik _  F-z 

2(z-u)F-u 


F-u 

F-z 


[21] 


The  closure  relations  for  <ya>,  <aV/x>  ,<Xaa>  and  <^aa*>  are 
derived  by  first  multiplying  eq.  20  by  X,,  £•£,  £a  and  £a*  respectively 
and  then  taking  the  ensemble  average.  Hypotheses  I,  II  and  III  are  used 
to  simplify  the  statistical  moments  of  the  form  <xU,£)  g(u,£)>  which 
are  operated  on  by  L(z,r;u,.s).  The  details  can  be  found  in  Refs.  18  and  20. 
We  simply  recall  the  results  here.  They  are 


<a V/Y>  =  -  K  <A> , 


[22] 


+There  is  an  error  in  the  definition  of  the  operator  L  in  Ref.  20.  The 
error  concerns  the  effect  of  the  curvature  F  and  disappears  in  the  limit 
of  collimated  beams.  Since  the  solutions  of  Ref.  20  were  calculated  in 
that  limit,  the  error  does  not  affect  the  results  and  conclusions. 
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<ya>  =  -  R  •  V  <A> , 


[23] 


<vaa>  =  -  1/2  P  •  V  <aa>  . 


[24] 


<vaa*>  =  -  1/2  Q  •  V  <aa*>  , 

rlz,  ~  » 


[25] 


where 


K  =  1/2  L(z,r;u,sJ  :  <v (z >r)xju ,£) >  ,  [26] 

R  =  L(z,r;  u,s)  <v(z  trjv(u,s]> ,  [27] 

P  =  &  ,  [28] 

rs*  ~ 

Q  =  Real  {£} .  [29] 


Equations  22-25  form  the  basis  of  our  model,  they  relate  the 
higher  order  unknown  moments  to  the  lower  order  moments  and,  thus,  close 
the  system  of  equations  for  <A>,  <aa>  and  <aa*>.  The  functional 
relationships  described  by  eqs.  22-25  are  consistent  with  the  physical 
interpretation  of  the  terms  they  model.  The  moment  <aV/v>  which  causes 
the  decay  of  the  coherent  amplitude  <A>  is  proportional  to  <A> ,  and, 
therefore,  it  survives  for  plane  and  spherical  waves  as  indeed  it  should. 
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But,  the  terms  <^a>,  <^aa>  and  <vaa*>,  responsible  for  beam  spreading, 
are  all  proportional  to  gradients  of  the  lower  order  mean-field  quanti¬ 
ties  and,  as  expected,  they  vanish  for  plane  and  spherical  waves.  It  is 
noteworthy  that  these  relations  are  analogous  with  expressions  success¬ 
fully  used  to  model  the  turbulent  momentum,  heat  and  mass  transport  terms 
in  turbulent  shear  flows  (see  for  example  Chap.  1,  5  and  6  of  Ref.  19). 
However,  they  were  not  simply  hypothesized  here,  they  were  derived  in  a 
fully  consistent  manner  from  the  governing  stochastic  wave  equation. 

Equations  26-29  show  that  the  proportionality  functions  K,  R,  P 

ft;  ~ 

and  are  all  integrals  of  the  covariance  of  the  angle-of-arrival 
<vi(z,r)vi(u,sj)>.  An  equation  for  v  is  obtained  by  subtracting  eq.  15 
from  eq.  3,  i.e. 


0 

—  v  +  <y>*Vv  +  v-V<V>  +  v*Vv  =  Vn/nD  • 
3z 


[30] 


Neglecting  the  nonlinear  stochastic  term  v/Vv,  which  is  valid  for  propa¬ 
gation  in  turbulence  since  the  fluctuating  ray  angle  |v|  remains  small 
throughout  the  propagation  range,  and  using  eq.  16  for  <V>,  we  have 


8  ~ 
—  v  +  - 

3 z  ~  (z-F) 


•Vv  + 


X, 

(7-~ 


V,n/n0. 


[31] 


Solving  along  the  average  geometrical  ray  that  passes  through  the  point 
(z ,t)  and  assuming  that  the  path  length  along  the  ray  can  be  approximated 
by  z  in  accordance  with  the  paraxial  approximation,  we  find 
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v(z,r) 


“/ 


dC 


u-n 


(z-F)  n0 


—  .  L'  =  —  l) 


F-z 


[32] 


Therefore,  the  covariance  of  v  is  given  by: 


cv(z,r)v(u,£)> 


(C-F)  q-F) 
(z-F)  (u-F) 


.  V 

~r 


<n 


r ' 


[33] 


Since  the  correlation  scale  of  the  index  covariance  function  is 
much  smaller  than  the  propagation  distances  of  interest,  eq.  33  can  be 
further  approximated  as  follows: 


±-[  dc-ii 

n|  J  (z- 


-F)‘ 


00 

-r'^r"  ^ 


F) (u-F) 


dr<n[c+T,  jr'  =  r"+A] 


134] 


The  vectors  r',  r"  and  A  are  defined  in  Fig.  1  which  illustrates  the 
geometrical  features  associated  with  the  definition  of  the  covariance 
function  of  v.  In  the  paraxial  approximation,  the  vector  £  is  given  by 


A 


F-A 

F-z 


r 


F-u 


[35] 


— "fom - t- 
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FIGURE  1  -  Configuration  of  the  rays  and  geometrical  definition  of  the 
symbols  for  the  calculation  of  the  covariance  function 


The  integral  over  x  has  been  worked  out  for  homogeneous  and  isotropic 
Kolmogorov  turbulence  (Refs.  1  and  2).  The  result  for  eq.  34  is 

f  ,  p,2  3.44  C2 

//  (z-F)  Cu-F)  ^  ‘  V*1  ■  S<<  ^  ' 


<v(z,jJ)CCu>s)> 


r  d?  - 

i  c 


/-  PI 2  1-62  C2 

77~~~  2  ai/3  [^a  +  A>>  to> 

z-F) (u-F)  no  A  2 


where  Cn  is  the  well-known  structure  constant  of  the  turbulent  refrac¬ 
tive  index  and  is  the  inner  scale  of  turbulence.  To  simplify  the 
algebra,  we  propose  the  following  simple  expression  connecting  the  two 
asymptotic  formula  of  eq.  36: 


2.43  C2  f 

^O.rjv/u,  s)>  =  - —  / 

n  l  J 


(z-F)  (u-F)  [A2  + 


j-TTT; 


35  £.)‘] 


nasBrawrasH 
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Before  substituting  eq.  37  for  the  covariance  of  ;v  in  eqs.  26  and  27 
for  the  proportionality  functions  K  and  _R,  it  is  essential  to  match  the 
geometry  associated  with  eq.  37  to  that  associated  with  the  definition 
of  the  operator  L.  This  is  best  explained  with  reference  to  Fig.  2. 

The  operator  L  was  derived  exactly  and  thus  applies  to  spherical  phase 
fronts.  In  particular,  the  argument  of  the  exponential  in  eq.  21  is 
the  phase  difference  between  wavelets  originating  at  A  and  B  and 
arriving  at  P,  i.e. 


A<f> 

1 


k  [CD  -  CB]  . 


[38] 


Since  all  angles  are  small  as  required  by  the  paraxial  approximation, 
we  have 


__2 

CD  ~  — -  ,  [39] 

2 (z-u) 


.2 


2 (F-u) 


[40] 


[41] 


Hence,  from  eqs.  38-40  it  follows  that 


A* 

1 


k 


F-z 


2(z-u) 


F-u 


F-z 


142] 
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which  is  identical  to  the  expression  in  eq.  21.  However,  as  shown  in 
Fig.  1,  the  geometry  used  to  derive  the  covariance  of  v_  is  such  that  the 
transverse  coordinate  is  not  taken  along  circular  arcs  centered  at  z  =  F 
but  rather  along  lines  normal  to  the  ray  passing  through  {z,rj.  Thus, 
the  argument  of  the  exponential  in  eq.  21  should  be  replaced  by  the 
phase  difference  between  the  wavelets  originating  at  A  and  C  respectively 
and  arriving  at  P,  plus  the  phase  difference  between  the  spherical  wave 
fronts  passing  through  points  A  and  C,  i.e. 


u 


FIGURE  2  -  Configuration  of  the  rays  and  geometrical  definition  of  the 
symbols  for  the  calculation  of  the  modified  operator 
Lfz,£;u,s) 
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Hence,  from  eqs.  39-41  and  43,  we  find 


k 

F  +  z  -2u 

F-u 

2  (z  -u) 

F-u 

F-z 

[44] 


Therefore,  eqs.  26  and  27  are  slightly  modified  as  follows: 


2_  oo 

K  =  ~  f  — —  if  d2s  V  7  :  <v  (z  ,r)v  (u ,  sj  > 

4tt  7  (u-z) 


j  ik  F+z-2u  ,?| 
.  exp  < -  -  ' 

( 2(z-u)  (F-U2  I 


z  00 

R=  it  -iH—  ff  d2s  <y  (z  ,£)v  (u  ,s)> 
2ir  J  (u-z)  77 


)  ik  F+z  -2u  ,2 

exp  <  -  -  A 

(  2  (z  -u)  (F-O2 


[45] 


[46] 


Equation  37  shows  that  the  covariance  of  v  is  function  of  A2 
alone  in  the  transverse  plane.  Since  the  integral  operator  in  eqs.  45  and  46 
is  also  dependent  on  A2  alone,  the  resulting  proportionality  functions 
K  and  R  depend  on  propagation  distance  z  only. 


I 

J 
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Equation  37  is  substituted  for  the  covariance  of  v_  in  eqs.  45  and  46. 
The  integration  over  the  transverse  plane  can  be  performed  analytically 
to  give 


K(z) 


f  du  f  d,  .ffrVl 

nl  J  J  (F-z)  (F+z-2u) 


[47] 


R(z) 


2.43  C2 

iTT7* 

1 


z  u 


0  o 


d? 


(F-O2  hCy) 
(F-z) (F+z-2u) 


[48] 


where 


G(y)  =  C-iy)2  -  3iY  ♦  [5/12  +  iY  +  3y2]— - FTT  r (5/6, -iY)| 

2  (-iY)Vb 


[49] 


-1Y 


H(y)  =  (-iY)  +  [5/6+iY]  - - - gyg-  T(5/6,  -iY)|»  [50] 


Y  = 


k£2 

.  1_ 


(F+z-2u)  (F-u) 


[51] 


2  (z-u)  (F-c)* 


and  where  t  =0.35  £0)6,is  the  unit  dyad  and  F(5/6,  -iy)  is  an  incomplete 
Gamma  function  as  defined  by  formula  6.5.20  of  Ref.  21.  The  functions 
G(y)  and  H(y)  are  plotted  in  Figs. 3  and  4  respectively. 
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The  closure  eqs.  22-25  together  witn  the  eqs .  47  and  48  for  the 
proportionality  functions  K(z)  and  R(zJ  are  sufficient  to  solve  for  the 
average  irradiance.  However,  to  calculate  the  irradiance  variance 
profiles,  we  need  a  relation  for  the  moment  <Va*Va>.  The  latter  was 
derived  in  a  semi -empirical  fashion  in  Ref.  20.  As  mentioned  in  the 
introduction,  a  more  theoretical  treatment  of  <£a-7,a>  would  require 
solving  for  the  covariance  function  of  the  amplitude,  hence  recourse  to 
mult idimensional  partial  differential  equations.  This  is  beyond  the 
scope  of  the  present  report  and,  therefore,  the  result  of  Ref.  20 
multiplied  by  a  geometrical  factor  to  account  for  focalization  will  be 
simply  taken  here,  i.e. 


<Va* Va> 


1  o 

—  Vz<aa> 
4 


c 

F-z 


C 

F-z 


-  <aa> ;  z  <  zA  , 


kz 


<aa>;  z  >  z. 
2  A 


[52] 


where  C  is  a  universal  empirical  constant  whose  best  rounded  value  is 
C=0.5-0.01  i^  and  z^  is  a  length  scale  defined  in  the  next  section. 
Finally,  by  virtue  of  the  approximation  of  quasi -homogeneity  and  quasi¬ 
isotropy,  i.e.  eq.  14b, and  from  standard  diffraction  results,  we  have 


t  In  Refs.  18  and  20,  the  constant  was  chosen  equal  to  0.5  -  0.01  i.  The 
effect  on  the  resulting  solutions  for  the  irradiance  standard  deviation 
is  less  than  10%.  These  differences  in  the  choice  of  C  are  explained 
by  the  different  methods  of  evaluating  the  functions  K(z)  and  R(z). 

In  Ref.  18,  we  limited  ourselves  to  real  and  purely  empirical  expres¬ 
sions;  in  Ref.  20,  we  used  asymptotic  formulas  only;  but,  here,  we 
calculate  K(z)  and  R(z)  from  the  full  integral  eqs.  47-48. 
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—  [ <a*V2a>  -  <aV2a*>]  ~  —  — £ —  V2< aa*>  .  [53] 

2k  2  (F-z) 


where  G  is  the  dif fractional  far-field  half-angle  beam  divergence.  For 
a  Gaussian  beam  with  an  irradiance  e-folding  radius  w  ,  0  =  1/kw^. 


4.0  COMPARISON  OF  SOLUTIONS  WITH  DATA 

After  substitution  of  eqs.  22-25,  52  and  53  for  the  higher  order 
moments  in  eqs.  10-12,  we  obtain  a  closed  set  of  partial  differential 
equations  for  <A>,  <aa>  and  <aa*>.  In  nondimensional  form,  these 
equations  are 


+  —  — Va>  -  V2  <A>-  b/?(n)V2<A>  =  - 


3n  n-f  3p  / 


-  -  £(n)<A>, 
(n-f)  2 


[54] 


3n  n-f  3p  / 


<aa>  -  —  v2  <aa>  -  —  /?(n)v2  <aa>  +  iT(n)<aa> 
4  p  2  p 


[55] 


2  <  cl  3.  ^ 

-  +  n)<A><A>  +  2b  /?(n)  V  <A>*V  <A>  , 

(n-f)  ~  ~D 


3  ,  P  3 

3n  n-f  3p 
2<aa*> 


)  <aa*>  -  ^  2  v- 

/  (f-n)  2  " 


V‘  <aa*>  -  —  Rea  1  [  (n)  ]  <aa*> 

2  ‘  D 


[56] 


+  Real [tf(n) ] <A><A>*  +  2b  Real[fi(n)lV  <A>*V  <A>* 
(n-f)  -f  "P 
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where  n  and  p  are  nondimensional  variables 


n  =  z/z^  and  p  =  r/w0  > 


[57a, b] 


and  b  and  f,  similarity  parameters 


b  =  2  /kw2  and  f=  F/z.  . 
a  o  A 


[58a,b] 


The  length  scales  w  and  z  are  respectively  the  representative  radius 
O  A 

of  the  unperturbed  laser  beam  and  the  fading  distance  of  the  average 
amplitude  given  by 


z  =  1/c12/ll  ,7/H 

A  n 


Finally,  the  proportionality  functions  are 


*n)  -  f  dy  /  dx  - 

no  J  J  (f-n)2  (f+n-2y) 


o  o 


n  y 


=  T7 \  f  ^  j  dx 


o  o 


(f-x)2H(y) 
(f-n) (f+n-2y) 


9 


[61] 
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T(n) 


(f-n/2)  i 
(f-n)  n 


; 

(f-n) 


n  <  1  , 


with 


no  (f+n-2y)  (f-y) 
2 (n-y)  (f-x)2 


which  gives  rise  to  a  third  similarity  parameter 


[62] 


[63] 


no  =  k  £2/ zA 


[64] 


However,  for  n>>no,  the  solutions  become  virtually  independent  of  no. 

Equations  54-56  form  a  closed  set  of  partial  differential 
equations,  in  three-dimensional  space  only,  for  the  solution  of  the 
statistical  moments  <A>,  <aa>  and  <aa*>.  General  solutions  to  these 
equations  can  be  worked  out  in  terms  of  Green's  functions  but  this  leads 
to  complicated  multiple  integrals  that  are  very  troublesome  and  lengthy 
to  evaluate.  We  find  it  much  more  convenient  to  solve  the  finite 
difference  version  of  eqs.  54-56.  This  approach  is  straightforward, 
involves  no  particular  difficulty  except  in  the  near  vicinity  of  n  =  f, 
and  is  applicable  to  arbitrary  beam  profiles. 

The  quantities  of  interest  are  the  average  irradiance  and  the 
irradiance  variance.  The  average  irradiance  is  defined  by 


,.u.  ,  "  SaiiMK*— *»-<■■■>  ■ 
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<1>  =  <A>  <A>*+<aa*>  ,  [65] 

and  the  irradiance  variance  by 

cr j  =  <(I-  <I>)2>  .  [66] 

In  terms  of  amplitude  moments,  we  have 

a2  =  <aa*aa*>  -  <aa*>2  +  2<A  ><aa*a*> 

+  2<A>*<aaa*>  +  2<A><A>*<aa*>  +  <A><A><aa>*  [67] 

+  <A>*<A>*<aa>  . 


Therefore,  the  irradiance  variance  depends  on  statistical  moments  of 
orders  three  and  four.  Instead  of  trying  to  derive  a  closed  set  of 
equations  for  these  higher  order  moments,  we  prefer,  as  a  first  and 
simpler  attempt  of  solution,  to  make  a  hypothesis  regarding  the  probabi¬ 
lity  distribution  of  the  complex-amplitude  fluctuations.  IVe  have  shown 
in  Refs.  22  and  23  that  the  hypothesis  of  normal  distribution  for  the 
complex  amplitude  afz.r)  constitutes  a  consistent  approximation  over 
the  complete  propagation  range.  The  assumption  of  normal  statistics 
for  a(z,r)  gives 


<a*a*a>  =  <aaa*>  =  0  , 


[68] 
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<aa*aa*>  =  2<aa*>2  +  <aa><aa>*  •  [69] 


Hence,  eq.  67  becomes 


a2  =  <aa*>2  +  <aa><aa>*  +  2<A><A>*<aa*> 

[70] 

+  <A><A><aa>*  +  <A>*<A>*<aa>  . 


Therefore,  both  the  average  irradiance  and  the  irradiance  variance  can 
be  calculated  from  the  solutions  of  eqs.  54-56.  It  is  noteworthy  that 
the  implications  of  normal  statistics  for  a(z,r)  are  explicitly  and 
completely  specified  by  eqs.  68  and  69.  Most  importantly,  this  hypothesis 
was  not  needed  to  derive  the  closure  relations  which  form  the  basis  of 
our  model. 

To  verify  our  model,  a  laboratory  experiment  was  designed  which 

allows  complete  control  over  the  turbulence  parameters  [Ref.  24).  In 

short,  the  refractive-index  turbulences  are  produced  by  creating  an 

unstable  vertical  temperature  gradient  in  a  tank  filled  with  water. 

This  is  simply  done  by  heating  the  water  at  the  bottom  of  the  tank  and 

cooling  it  at  the  top.  The  tank  is  1.5  m  long,  0.6  m  deep  and  0.4  m 

wide.  The  propagation  path  is  increased  by  folding  the  beam  lengthwise. 

Typically,  the  index  structure  constant  C  is  10~4m"1//3. 

n 

Figures  5  and  6  show  the  measured  normalized  average  beam  radius 
w/wo  plotted  versus  propagation  distance  for  collimated  (f  =  <*>)  Gaussian 
beams  with  w0  equal  to  6.0  mm  (b=0 . 0016)  and  3.0  mm  (b=0.0064)  respecti¬ 
vely.  The  data  were  taken  by  slowly  traversing  a  point  detector  across 
the  beam  through  its  center.  The  radius  w  was  calculated  from  the 
measured  average  irradiance  profile  <I(r)>  according  to 


a 
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FIGURE  5  -  Normalized  beam  radius  plotted  versus  normalized  propagation 
distance;  b  =  0.0016,  f  =  °°,  no  =  7.0.  □  :  data  measured 

in  laboratory  turbulence.  -  :  finite-difference 

solution  of  eqs.  54  and  56. 


FIGURE  6  -  Normalized  beam  radius  plotted  versus  normalized  propagation 
distance;  b  =  0.0064,  f  =  no  5  7.0.  □  :  data  measured 

in  laboratory  turbulence.  - :  finite-difference 

solution  of  eqs.  54  and  56. 
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w 


2 


<1 (r) >dr 


o 


r  < I  (r)>  dr . 


[71] 


The  solid  curves  represent  the  theoretical  solutions  calculated  from 
eqs .  54  and  56.  The  predicted  beam  spread  is  smaller  than  the  measured 
spread  by  about  10“.  However,  the  experimental  calculation  of  the  beam 
radius  becomes  quite  sensitive,  with  increasing  radius,  to  errors  in 
the  measured  irradiance  as  indicated  by  the  presence  of  r  to  the  third 
power  in  the  integral  of  the  numerator  of  eq.  71.  Since  irradiance 
errors  at  low  level,  such  as  in  the  wings  of  the  beam,  are  generally 
toward  higher  values,  and  overestimation  of  w2  appears  likely.  There¬ 
fore,  we  conclude  that  the  turbulence-induced  beam  spread  is  satisfac¬ 
torily  predicted  by  our  model.  In  particular,  the  influence  of  the 
similarity  parameter  b  =  z^/kws  is  well  verified. 

Figures  7  and  8  show  similar  results  for  beams  focused  at 
f  =  3.13  and  f  =  4.41  respectively.  The  data  were  measured  from 
exposed  photographs  of  the  beam.  The  images  were  digitized  and  the 
radii  calculated  with  the  formula 


w 


2 


dxdy (xz 


'2) 


<1  (x,y)> 


dxdy  <1 (x ,y) >  . 


[72] 


In  this  case,  the  agreement  between  the  data  and  the  predictions  of  our 
model  is  excellent.  The  average  beam  radius  passes  through  a  minimum 
which  is  considerably  greater  than  the  diffraction  limit  value  equal  to 
(bf)  and  this  occurs  at  a  propagation  distance  appreciably  shorter  than 
the  focal  length  f.  The  magnitude  of  this  effect  and  its  dependence  on 
the  parameter  f  are  well  predicted  by  our  model  as  illustrated  in 
Figs.  7-8. 
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FIGURE  7 


FIGURE  8 


-  Normalized  beam  radius  plotted  versus  normalized  propagation 
distance;  f  r  3.13,  b  =  0.0023,  no  =  8.5.  □  :  data  measured 

in  laboratory  turbulence.  - :  finite-difference 

solution  of  eqs.  54  and  56. 


-  Normalized  beam  radius  plotted  versus  normalized  propagation 
distance;  f  =  4.41,  b  =  0.0021,  no  =  8.5.  □  :  data  measured 

in  laboratory  turbulence.  - :  finite-difference 

solution  of  eqs.  54  and  56. 
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FIGURE  9  -  Normalized  irradiance  standard  deviation  plotted  versus 

normalized  propagation  distance.  Data  measured  in  laboratory 
turbulence  on  the  axis  of  a  25-mm  diameter  collimated  (f  =  °°) 
beam.  0:rio  =  4.0,  b  =  0.0020;  A  :  n0  =  5.3,  b  =  0.0015; 

V  :  no  =  7.0,  b  =  0.0011.  - :  finite-difference  solution 

of  eqs .  54  -  56,  curve  No  1  for  n<>  =  4.0,  b  =  0.0020  and 
curve  No  2  for  n<>  =  7.0,  b  =  0.0011. 


FIGURE  10  -  Normalized  irradiance  standard  deviation  versus  normalized 
propagation  distance.  Data  measured  in  the  atmosphere  and 
reprinted  from  Fig.  19  of  Ref.  7.  - :  finite-differ¬ 

ence  solution  of  eqs.  54-56  for  plane  waves  (f  =  °°,  b  =  0) . 
- :  asymptotic  solutions  derived  in  Ref.  7. 
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Figure  9  compares  the  irradiance  standard-deviation  data  measured 
on  the  axis  of  collimated  beams  with  the  prediction  of  our  model.  These 
results  served  to  determine  the  empirical  constant  C  of  eqs.  52  and  62.  To 
verify  that  C  is  universal,  the  atmospheric  data  of  Fig.  19  of  Ref.  7  are 
plotted  in  Fig.  10  and  compared  with  our  plane-wave  solution.  The 
agreement  is  excellent.  Therefore,  the  results  of  Figs  9-10  show  that 
our  model  works  very  well  under  experimental  conditions  that  are  as  far 
apart  as  can  be  expected  in  practice.  Indeed,  the  propagation  distances 
and  the  turbulence-strength  parameters  C^'s  differ  by  more  than  two 
orders  of  magnitude.  Therefore,  the  empirical  constant  C  can  in  fact 
be  considered  universal  with  good  accuracy. 

For  comparison  purposes,  we  have  shown  in  Fig.  10  the  asymptotic 
solutions  derived  in  Ref.  7  which  are  representative  of  the  present 
state  of  the  art.  We  therefore  see  that  our  solutions  agree  much  more 
closely  with  data,  particularly  with  regard  to  the  slow  return  of  the 
data  to  the  asymptotic  saturation  level  of  unity. 


5.0  SIMPLIFIED  SOLUTIONS 

Although  eqs.  54-56  together  with  the  constitutive  eqs.  60-62 
are  easily  solved  by  numerical  techniques,  it  is  useful  to  consider 
situations  where  approximate  analytic  solutions  can  be  worked  out. 
Despite  the  necessary  restrictive  conditions  under  which  these  solutions 
apply,  the  analytic  expressions  give  helpful  information  regarding  the 
influence  of  the  various  parameters  and  the  magnitude  of  the  turbulence- 
induced  phenomena.  One  such  simplification  is  afforded  by  the  study  of 
untruncated  Gaussian  beams. 

5. 1  Average  irradiance 


For  an  initially  Gaussian  beam  with  the  normalized  average 
amplitude  given  by 
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<A> 


n=o 


=  e 


-p2/: 


[73] 


it  is  straightforward  to  show  that  the  solution  of  eq.  54  is 


<A>  =  g(n) 


e-o2/2o2  (n) 


[74] 


where 


g(n) 


[75] 


a2(n) 


(f-n)2 

f2 


bi  nincil  +  2b  f 

£  Jo 


---  . 

(f-O2 


[76] 


If  turbulence- induced  beam  spreading  dominates  over  diffraction-spreading 
and  if  the  condition 


[ Imag(a2) ] 2/[Real (a2) ] 2  << 


[77] 


is  satisfied,  which  occurs  in  most  applications,  eq.  56  for  <aa*>  can 
be  solved  analytically  and  we  find  for  the  average  irradiance 
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<I>  =  {1  -  h(n)[l  -  P2/w“(n)]}  e 


p2/w2 (n) 


[78] 


w2  (n) 


where 


h(n) 


2b  (f-n) ' 

w2(n)  f2  JQ  wz(0 


n  c 

/"  — -  Real  [/? CO ]  exp {-  f  Real  [xfC)  ] d£ ) , 

/  w2CO  7 


[79] 


w2(n)  =  +  2hJ  Real  [7?(?)]d?. 


f2 


[80] 


lf-0' 


Equation  78  shows  that  the  turbulence-spread  Gaussian  beam  does 
not  truly  conserve  its  Gaussian  shape;  too  much  irradiance  is  scattered 
in  the  wings  as  inferred  by  the  p2-term  multiplying  the  exponential 
function.  However,  in  many  practical  applications ,  the  similarity 
parameter  b  is  much  smaller  than  unity  and  the  function  h(n)  can  be 
neglected  in  eq.  78  to  give 

-p2/w2(n) 

<I>  =  £ -  .  [81] 

w2(n) 


An  important  simplification  can  be  made  to  the  proportionality 

function  I?(n)  if  we  assumed  that  no  «  1.  This  condition  is  almost 

always  satisfied  for  atmospheric  propagation  of  optical  and  infrared 

beams.  For  example,  at  X  =  0.63  pm,  no  is  equal  to  about  10~3  for  very 

strong  turbulence  (C  =  iq-6  -l/3  )  and  smaller  for  weaker  turbulence. 

n 

Under  this  condition,  the  asymptotic  formula 


4 
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l 


lim  H(y)  =  -  r(5/6)(-iy)'1/6  f82l 

y-*-0  6 


can  be  substituted  for  H(y)  in  eq.  61  and  we  find 


Real[ft(n)] 


0.738 


I 


Cf8/3-(f-y)a/31(f-y) 


dy 


8/3-i 


,i/6 


(f-n)(f+n 


,  ,5/6,  ,1/6 

-2y)  (n-y) 


[83] 


This  integral  cannot  be  worked  out  in  terms  of  elementary  transcendental 
functions.  However,  for  the  following  discussion,  it  suffices  to  consider 
collimated  beams,  i.e.  f  =  “>,  which  yields 


Real ]  =  1.288  r)11/6 


[84] 


Upon  substituting  eq.  84  for  i?(n)  in  eq.  80,  we  therefore  have  for 
collimated  beams 


[85] 

Thus,  the  turbulence-induced  beam  spread  predicted  by  our  model  and 
valid  under  the  conditions  [77]  and  q0  <<  1  is  given,  in  dimensional 
variables,  by 


w2 (n)  =  1  +  0.909  b  n17/6 
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p2(z)  =  w2  +  0.909  k 


1/6 


C2 

n 


,17/6 


[86] 


The  beam-spread  expression  most  often  used  in  the  literature 
is  that  derived  by  Brown  (Ref.  25)  and  Yura  (Ref.  26) .  It  is  based  on 
the  reasoning  that  the  turbulence-induced  beam  spread  can  be  no  smaller 
than  that  produced  by  a  transmitter  with  an  aperture  equal  to  the 
oherence  diameter  p0  =  (0.545  k2zC2/n2)  The  formula  given  by 

Yura  is 


2  2  1  07  1.2/5  -16/5 

Py  =  o  +  i ' 93  k  z 


(Cn/n,) 


12/5 


[87] 


We  note  that  both  ours  and  Brown's  or  Yura's  expressions  depend 
on  the  same  three  parameters  k,  C^/no  and  z.  In  particular,  they  are 
independent  of  beam  or  aperture  size  which  confirms  the  experimental 
observation  that  turbulence  spreading  cannot  be  corrected  by  increasing 
the  transmitter  aperture.  However,  the  functional  dependence  on  k, 

C^/no  and  z  is  slightly  different  in  both  cases.  Our  expression  predicts 
a  weaker  dependence  on  all  three  parameters. 

Systematic  turbulence  beam-spread  measurements  were  performed 
by  Dowling  and  Livingston  (Ref.  27).  For  total  far-field  turbulence- 
induced  beam  divergence,  they  proposed  the  following  empirical  regression 
formula 


02dl  =  2.9  x  10'6  k1/3  Z  (Cn/n0)6/5. 


[88] 
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Table  l  compares  the  predictions  obtained  from  eqs.  86  and  87,  i.e. 


’.27  k1/6  Z5/6  (C  /nD)2  , 


[89] 


02 

Y 


15.44  k2/S  Z6/5  (Cn/no)12/S 


[90] 


with  those  of  eq.  88  for  the  upper,  middle  and  lower  values  reported 
in  Ref.  27.  Our  model  gives  a  somewhat  more  consistent  fit  to  the 
empirical  formula.  The  discrepancies  come  mainly  from  our  stronger 
dependence  compared  with  that  of  eq.  88.  However,  it  is  shown  in  the 
same  Ref.  27  that  the  optical  data  (X  =  0.63  ym)  are  also  well  fitted 
by  a  linear  regression  against  C2z  which  has  the  same  C  power  as  our 
model  and  only  a  difference  of  1/6  in  the  power  of  propagation  distance. 
Therefore,  the  beam-spread  formula  derived  from  the  propagation  model 
of  this  report  appears  to  be  more  precise  than  the  generally  used 
formula  of  Brown  and  Yura,  although  more  data  with  greater  parameter 
ranges  are  still  needed. 
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TABLE  I 

Comparison  of  the  theoretical  beam-spread  formula  of  this  report 
(0^,  eq.  89)  and  that  of  Ref.  26  (Oy ,  eq.  90)  with  the  empirical 
regression  formula  (0^,  eq-  88)  proposed  in  Ref.  27.  The  propagation 
distance  z  is  1750  m  and  the  C^-values  correspond  to  the  upper,  middle 
and  lower  turbulence  levels  reported  in  Ref.  27. 


A  =  0.63  urn 

■E 

10.6  ym 

c 

n 

®DL 

0B 

eY 

SDL 

0B 

0Y 

,  -1/3-. 

(m  ) 

(yrad) 

(yrad) 

(yrad) 

(yrad) 

(yrad) 

(yrad) 

8  X  10'7 

229 

185 

420 

143 

147 

239 

5  X  10"7 

173 

116 

239 

108 

92 

136 

1  x  io"7 

66 

23 

35 

41 

18 

20 
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5.2  Irradiance  variance 


Under  the  same  simplifying  assumptions  as  used  to  derive  eq.  81 
for  the  average  irradiance,  it  is  easy  to  show  from  eqs.  54-56  and  70  that 
the  corresponding  irradiance  variance  profile  has  the  form 


,2  _ 


=  F (n)  e 


-p2/w2(n) 


[91] 


where  F(n)  contains  integrals  that  cannot  be  reduced  to  elementary 
transcendental  functions.  However,  the  important  conclusion  from 
eqs.  91  and  81  is  that  the  normalized  irradiance  variance,  o2/<I>2,  is 
independent  of  lateral  position  across  the  beam  but  varies  with  propa¬ 
gation  distance  only.  Within  the  experimental  scatter,  this  is  well 
verified  in  our  simulated  turbulent  medium. 

For  plane  waves,  the  function  F(n)  of  eq.  91  reduces  to  a 
simple  asymptotic  formula  valid  for  n  <  1.  The  solution  for  the  irra¬ 
diance  variance  is  thus 


<I>2 


=  2.76 


k7/6  Z1176  cVni 
n 


[92] 


Except  for  the  numerical  constant,  this  expression  is  identical  with  the 
perturbation  result  given  by  Tatarskii  (Ref.  2),  i.e. 


1.23  k7/6  Z11/6  C2/n2 
n 


[93] 
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Equation  92  corresponds  to  the  weak-scintillation  asymptote  of  the  plane- 
wave  solution  shown  as  the  solid  curve  in  Fig.  10  while  Tatarskii's 
solution,  eq.  93,  is  plotted  as  the  broken  line  through  the  origin  in 
the  same  figure.  We  thus  see  that  either  solution  is  equally  justifiable 
on  the  basis  of  the  low-scintillation  data  of  Fig.  10.  Therefore,  we 
conclude  that  the  propagation  model  proposed  in  this  report  is  in  good 
agreement  with  the  results  of  the  classical  perturbation  theory. 


6.0  CONCLUSION 

The  mathematical  model  of  laser-beam  propagation  in  turbulence 
originally  developed  in  Refs.  18  and  20  was  extended  to  include  focal- 
ization.  The  model  is  in  the  form  of  a  closed  set  of  three  simulta¬ 
neous  second-order  partial  differential  equations  for  the  complex- 
amplitude  moments  of  order  one  and  two.  The  important  and  novel 
characteristic  of  these  equations  is  that  they  are  numerically  tractable 
and  uniformly  valid  for  arbitrary  scintillation  levels. 

The  solutions  were  calculated  with  a  finite-difference  algorithm. 
The  agreement  with  measured  average  irradiance  profiles  for  collimated 
and  focused  beams  is  excellent.  The  predicted  irradiance  variance 
clearly  shows  super saturation  and  is  well  verified  by  data  ranging  from 
the  perturbation  to  the  asymptotic  saturation  regimes.  An  analytic 
expression  for  the  turbulence -induced  spreading  of  Gaussian  beams  was 
derived.  Comparison  with  data  showed  that  it  is  somewhat  more  precise 
than  the  commonly  used  formula  based  on  the  concept  of  the  coherence 
diameter  p0.  Finally,  the  irradiance  variance  solution  in  the  weak- 
scintillation  limit  agrees  with  the  classical  results  of  the  perturbation 
theory. 
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Extension  of  the  model  to  media  with  varying  and/or  nonlinear 
average  refractive  index  is  straightforward.  In  particular,  the  thermal - 
blooming  effect  in  the  presence  of  turbulence  could  be  calculated  with 
only  a  small  increase  in  computation  time  compared  with  the  non-turbulent 
case.  The  model  could  also  be  used  to  simulate  the  propagation  of  beams 
corrected  by  adaptive  systems.  The  principal  requirement  would  be  to 
rederive  the  solutions  for  the  functions  K(z)  and  R( z)  with  the  proper 
boundary  conditions  on  the  covariance  of  the  angle-of-arrival . 
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